#ifndef UTILITY_H
#define UTILITY_H

#include "CUDA_vec.h"
#include <stdio.h>

template<typename T> T dotProd(const CUDA_vec<T> & a, const CUDA_vec<T> & b){ return a.dotProd(b); }

template<typename T> T norm(const CUDA_vec<T> & a) { return (T)sqrt(a.dotProd(a)); }

void print_vector(const CUDA_vec<float> & a);

void print_2Dvector(const CUDA_vec<float> &a);

void print_CUDA_err_msg(const char * msg);

int max_threads_per_block(int devnum);

dim3 max_dim_block(int devnum);

dim3 max_dim_grid(int devnum);

dim3 optimal_block_size(int devnum, int dimension);

void print_dim(dim3 d);

template<typename T>
void print_device_array(const T* device_array, int size);

template<typename Matrix>
void compute_residual(const Matrix & A,
		const CUDA_vec<float> & x,
		const CUDA_vec<float> & b,
		CUDA_vec<float> & r )
{
	A(x, r);
	r.AXPBY(-1, 1, b);
}

template<typename Matrix>
float residual_norm(const Matrix & A, const CUDA_vec<float> & x, const CUDA_vec<float> & b)
{
	CUDA_vec<float> r(b.getSize());
	compute_residual(A,x,b,r);

	return norm(r);
}


#endif
